
#define ALIGN     16
#define ROUND(n)  (((size_t)(n)+ALIGN-1)&(~(size_t)(ALIGN-1)))

// This kernel is executed with a very high number of threads = total number of vertices * 3
__kernel
void eating(__global float *vbo, float dy, unsigned vertices_per_atom)
{
   int index = get_global_id(0);
   int tpa = 3 * vertices_per_atom; // threads per atom
   // ghost or pacman?
   float enable = (index / tpa) % 2 ? 0.0f : 1.0f;
   // y coordinate?
   float is_y = ((index % 3) == 1) ? 1.0f : 0.0f;
   // vertex belongs to upper part of sphere?
   float is_up = (index % tpa >= tpa/2) ? 1.0f : -1.0f;
   
   vbo[index] += dy * is_y * is_up * enable;
}

__kernel
void move_vertices(__global float *vbo, __global float *speed, unsigned vertices_per_atom)
{
  int index = get_global_id(0);
  int tpa = 3 * vertices_per_atom;
  int numero_atom = index/tpa;

  vbo[index] += speed[(index % 3) * ROUND(get_global_size(0)/tpa) + index/tpa];
}

__kernel
void update_position(__global float *pos, __global float *speed, unsigned N)
{
  int index  = get_global_id(0);
  pos[index] += speed[index];
  pos[index + ROUND(N)] += speed[index + ROUND(N)];
  pos[index + 2*ROUND(N)] += speed[index + 2* ROUND(N) ];
}

__kernel
void border_collision(__global float *pos, __global float *speed, __constant float *min, __constant float *max, float radius, unsigned N)
{ 
  int index = get_global_id(0);
  if(pos[index] - radius <= min[0])
    speed[index] = -speed[index];
  else if (pos[index] + radius >= max[0])
    speed[index] = -speed[index];

  if(pos[index + ROUND(N)] - radius <= min[1])
    speed[index + ROUND(N)] = -speed[index + ROUND(N)];
  else if (pos[index + ROUND(N)] + radius >= max[1])
    speed[index + ROUND(N)] = -speed[index + ROUND(N)];

    if(pos[index + 2*ROUND(N)] - radius <= min[2])
    speed[index + 2*ROUND(N)] = -speed[index + 2*ROUND(N)];
  else if (pos[index + 2*ROUND(N)] + radius >= max[2])
    speed[index + 2*ROUND(N)] = -speed[index + 2*ROUND(N)];


}

void atom_colision_line(int index,float * pos, float * speed, float radius, int N){
	
	float3 myPos;
	myPos.x = pos[index];
	myPos.y = pos[index + ROUND(N)];
	myPos.z = pos[index + 2 * ROUND(N)];
	
	float3 otherPos;
    for (int i = 0; i < index; i++) {
		otherPos.x = pos[i];
		otherPos.y = pos[i + ROUND(N)];
		otherPos.z = pos[i + 2 * ROUND(N)];
		if(distance(myPos, otherPos) <= 2 * radius) {
			speed[index] = 0;
			speed[index + ROUND(N)] = 0;
			speed[index + ROUND(N) * 2] = 0;
			speed[i] = 0;
			speed[i + ROUND(N)] = 0;
			speed[i+ 2* ROUND(N)] = 0;
		}
	}
	
}

// V1
//~ __kernel
//~ void atom_collision(__global float *pos, __global float *speed, float radius)
 //~ {
  //~ int index = get_global_id(0);
  //~ int N = ROUND(get_global_size(0));
  //~ float3 myPos;
  //~ myPos.x = pos[index];
  //~ myPos.y = pos[index + ROUND(N)];
  //~ myPos.z = pos[index + 2 * ROUND(N)];
  //~ float3 otherPos;
  //~ for (int i = 0; i < index; i++) {
    //~ otherPos.x = pos[i];
    //~ otherPos.y = pos[i + ROUND(N)];
    //~ otherPos.z = pos[i + 2 * ROUND(N)];
    //~ if(distance(myPos, otherPos) <= 2 * radius) {
      //~ speed[index] = 0;
      //~ speed[index + ROUND(N)] = 0;
      //~ speed[index + ROUND(N) * 2] = 0;
      //~ speed[i] = 0;
      //~ speed[i + ROUND(N)] = 0;
      //~ speed[i+ 2* ROUND(N)] = 0;
    //~ }
  //~ }
//~ }

//#######################  V2  ############################
//~ __kernel
//~ void atom_collision(__global float *pos, __global float *speed, float radius)
//~ {
	//~ int index = get_global_id(0);
    //~ int N = ROUND(2*get_global_size(0));
	//~ 
	//~ float3 myPos;
	//~ myPos.x = pos[index];
	//~ myPos.y = pos[index + ROUND(N)];
	//~ myPos.z = pos[index + 2 * ROUND(N)];
	//~ 
	//~ float3 otherPos;
    //~ for (int i = 0; i < index; i++) {
		//~ otherPos.x = pos[i];
		//~ otherPos.y = pos[i + ROUND(N)];
		//~ otherPos.z = pos[i + 2 * ROUND(N)];
		//~ if(distance(myPos, otherPos) <= 2 * radius) {
			//~ speed[index] = 0;
			//~ speed[index + ROUND(N)] = 0;
			//~ speed[index + ROUND(N) * 2] = 0;
			//~ speed[i] = 0;
			//~ speed[i + ROUND(N)] = 0;
			//~ speed[i+ 2* ROUND(N)] = 0;
		//~ 
		//~ }
	//~ }
	//~ 
	//~ index = N-index;
	//~ 
	//~ myPos.x = pos[index];
	//~ myPos.y = pos[index + ROUND(N)];
	//~ myPos.z = pos[index + 2 * ROUND(N)];
	//~ 
    //~ for (int i = 0; i < index; i++) {
		//~ otherPos.x = pos[i];
		//~ otherPos.y = pos[i + ROUND(N)];
		//~ otherPos.z = pos[i + 2 * ROUND(N)];
		//~ if(distance(myPos, otherPos) <= 2 * radius) {
			//~ speed[index] = 0;
			//~ speed[index + ROUND(N)] = 0;
			//~ speed[index + ROUND(N) * 2] = 0;
			//~ speed[i] = 0;
			//~ speed[i + ROUND(N)] = 0;
			//~ speed[i+ 2* ROUND(N)] = 0;
		//~ }
	//~ }
//~ }


// ############################  V3  #####################
__kernel
void atom_collision(__global float *pos, __global float *speed, float radius, int natoms)
{
	int index = get_global_id(0);
	int N = ROUND(natoms);
	
	float3 myPos;
	myPos.x = pos[index];
	myPos.y = pos[index + ROUND(N)];
	myPos.z = pos[index + 2 * ROUND(N)];
	
	for(int i = 0; i < index; i++){
		float3 otherPos;
		otherPos.x = pos[i];
		otherPos.y = pos[i + ROUND(N)];
		otherPos.z = pos[i + 2 * ROUND(N)];
		
		if(distance(myPos, otherPos) <= 2 * radius){
			
			float3 vect_i = normalize(otherPos-myPos);
			
			float3 mySpeed, otherSpeed;
			mySpeed.x = speed[index];
			mySpeed.y = speed[index + ROUND(N)];
			mySpeed.z = speed[index + 2 * ROUND(N)];
			
			otherSpeed.x = speed[i];
			otherSpeed.y = speed[i + ROUND(N)];
			otherSpeed.z = speed[i + 2 * ROUND(N)];
			
			float3 vect_j;
			vect_j.x = 0;
			vect_j.y = -vect_i.z;
			vect_j.z =  vect_i.y;
			
			vect_j = normalize(vect_j);
			
			float3 vect_k;
			vect_k.x = vect_i.y*vect_j.z - vect_i.z*vect_j.y;
			vect_k.y = vect_i.x*vect_j.z - vect_i.z*vect_j.x;
			vect_k.z = vect_i.x*vect_j.y - vect_i.y*vect_j.x;
			
			float tmp = 1 + vect_i.x;
			
			if(tmp < 0.1)
				tmp = 0.1;
			
			float3 Ml1;
			Ml1.x = vect_i.x;
			Ml1.y = vect_i.y;
			Ml1.z = vect_i.z;
			
			float3 Ml2;
			Ml2.x = -vect_i.y;
			Ml2.y =  vect_i.x + (vect_i.z*vect_i.z)/tmp;
			Ml2.z = -(vect_i.y*vect_i.z)/tmp;
			
			float3 Ml3;
			Ml3.x = -vect_i.z;
			Ml3.y = -(vect_i.y*vect_i.z)/tmp;
			Ml3.y = vect_i.x + (vect_i.y*vect_i.y)/tmp;
			
			float3 Ml1_t;
			Ml1_t.x =  vect_i.x;
			Ml1_t.y = -vect_i.y;
			Ml1_t.z = -vect_i.z;
			
			float3 Ml2_t;
			Ml2_t.x = vect_i.y;
			Ml2_t.y = vect_i.x + (vect_i.z*vect_i.z)/tmp;
			Ml2_t.z = -(vect_i.y*vect_i.z)/tmp;
			
			float3 Ml3_t;
			Ml3_t.x = 	vect_i.z;
			Ml3_t.y = -(vect_i.y*vect_i.z)/tmp;
			Ml3_t.z = 	vect_i.x + (vect_i.y*vect_i.y)/tmp;
			
			float3 tmp1, tmp2;
			tmp1.x = dot(Ml1, mySpeed);
			tmp1.y = dot(Ml2, mySpeed);
			tmp1.z = dot(Ml3, mySpeed);
			
			tmp2.x = dot(Ml1, otherSpeed);
			tmp2.y = dot(Ml2, otherSpeed);
			tmp2.z = dot(Ml3, otherSpeed);
			
			tmp 	= tmp1.x;
			tmp1.x	= tmp2.x;
			tmp2.x 	= tmp;
			
			speed[index] 				= dot(Ml1_t, tmp1);
			speed[index + ROUND(N)] 	= dot(Ml2_t, tmp1);
			speed[index + ROUND(N) * 2] = dot(Ml3_t, tmp1);
			
			speed[i] 				= dot(Ml1_t, tmp2);
			speed[i + ROUND(N)] 	= dot(Ml2_t, tmp2);
			speed[i + ROUND(N) * 2] = dot(Ml3_t, tmp2);
			
			pos[index] -= vect_i.x/100;
			pos[index + ROUND(N)] 	-= vect_i.y/100;
			pos[index + ROUND(N) * 2] -= vect_i.z/100;
			
			pos[i] 				+= vect_i.x/100;
			pos[i + ROUND(N)] 	+= vect_i.y/100;
			pos[i + ROUND(N) * 2] += vect_i.z/100;
		}
	}
			
}


__kernel
void gravity(__global float *pos, __global float *speed, float g)
{
	int index = get_global_id(0);
	int N = get_global_size(0);
	
	speed[ROUND(N) + index] -= g;
	
	
}


__kernel
void lennard_jones(__global float *pos, __global float *speed, float radius, int natoms)
{
	int index = get_global_id(0);
    int N = ROUND(get_global_size(0));

	float3 myPos;
	myPos.x = pos[index];
	myPos.y = pos[index + ROUND(N)];
	myPos.z = pos[index + 2 * ROUND(N)];
	
	float3 otherPos;
	float3 speed_tmp1 = {0, 0, 0};
	
	if(index < natoms){
		
	    for (int i = 0; i < index; i++) {
			otherPos.x = pos[i];
			otherPos.y = pos[i + ROUND(N)];
			otherPos.z = pos[i + 2 * ROUND(N)];
			
			float d = distance(myPos, otherPos);
			
			float3 vect_i = -normalize(otherPos - myPos);
					
			float eps = 0.005;
			//~ float sigma = 2.2 * radius;
			float sigma = 1.1;
			
			float3 F = { 0., 0., 0.};
			float3 x = { 1., 0., 0.};
			float3 y = { 0., 1., 0.};
			float3 z = { 0., 0., 1.};
			
			float tmp1 = pow(d,12.);
			float tmp2 = pow(d,6.);
			
			if(tmp1 < 0.001){
				tmp1 = 0.001;
				tmp2 = pow(1/0.001, 6.);
			}
			
			float f = 4 * eps * ( pow(sigma, 12.)/tmp1 - pow(sigma, 6.)/tmp2);
			F.x = f * vect_i.x;
			F.y = f * vect_i.y; 
			F.z = f * vect_i.z;
			
			barrier(CLK_GLOBAL_MEM_FENCE);
			
			speed[index] 				+= F.x;
			speed[index + ROUND(N)] 	+= F.y;
			speed[index + ROUND(N) * 2] += F.z;
			
			speed[i] 				-= F.x;
			speed[i + ROUND(N)] 	-= F.y;
			speed[i + ROUND(N) * 2] -= F.z;
			
		}
	
		barrier(CLK_GLOBAL_MEM_FENCE);

	}
}
